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We present a numerical method to compute the Landauer conductance of non-interacting two- 
dimensional massless Dirac fermions in disordered systems. The method allows for the introduction 
of boundary conditions at the ribbon edges and to account for an external magnetic field. By 
construction, the proposed discretization scheme avoids the fermion doubling problem. The method 
does not rely on an atomistic basis and is particularly useful to deal with long-range disorder, whose 
correlation length largely exceeds the underlying material crystal lattice spacing. As an application, 
we study the case of monolayer graphene sheets with zigzag edges subjected long range disorder, 
which can be modeled by a single-cone Dirac equation. 
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I. INTRODUCTION 

The growing research activity on the physical proper- 
ties of graphen^ as well as on topological insulatorPE! 
has increased the interest in methods to solve the two 
dimensional massless Dirac equation for various settings. 
In pristine graphene, two Dirac cones placed at the cor- 
ners of the Brillouin zone, frequently called valleys, pro- 
vide a good description of the electronic single-particle 
low-energy spectrum. For topological insulators, two- 
dimensional massless Dirac excitations appear as edge 
states between a non trivial tridimensional (3D) topolog- 
ical insulator and a trivial vacuum. 

In distinction to analytical studies, that describe the 
single-particle properties of graphene by an effective two- 
dimensional Dirac Hamiltonian, almost all numerical 
analysis employ an atomistic basis from the onset (for 
reviews, see for instance Refs. 14151) . Numerical studies of 
the electronic transport in graphene sheets and ribbons 
typically use the tight-binding approximation to express 
the relevant Green's functions in terms of atomic site 
positions.'SHU This basis choice has the advantage of sim- 
plicity, but its computational cost increases very fast with 
system size. For the cases where short range disorder is 
important and an accurate description of the system at 
small length scales is necessary, the use of an atomic basis 
is not only the simplest, but also the most natural choice. 
The picture is very different in systems characterized by 
long range disorder and/or subjected to a smooth back- 
ground electrostatic p oten tial. In these cases, interval- 
ley scattering is absenlPSl and the system description in 
terms of a single Dirac cone is very accurate. Here, a nu- 
merical method that starts with the Dirac Hamiltonian is 
certainly computationally more efficient and very likely 
to be more insightful. The development of such method 
is the main goal of this paper. 

Topological insulator J2EI have rekindled the interest on 
the research related to concepts of topological order, of- 



ten used to characterize the singular behavior of the in- 
teger and fractional quantum Hall effects. Topological 
insulators are mainly characterized by the presence of 
chiral boundary states. It was recently discovered a se- 
ries of next generation 3D topological insulators. These 
materials exhibit a very simple band structure with a 
single Dirac cone. The non trivial topological nature of 
these materials makes their boundary states more ro- 
bust against disorder than the corresponding ones in 
graphene. In lattice quantum field theory the introduc- 
tion of boundary conditions has been guided by the in- 
dex theorem^ and by the phenomenological presence of 
Adler-Bell-Jackiw (chiral) anomalied^ observed in high 
energy experiments. In this context, it has been shown 
that nonlocal boundary conditions are necessary to re- 
produce the observed anomaly. For graphene and topo- 
logical insulator s the boundary conditions are local and 
well understood.'^^^ 

In this study we combine two finite difference 
methods ,^^^1111 developed in the context of lattice gauge 
theory, to numerically compute the scattering matrix of 
two-dimensional massless Dirac particles. Our method 
avoids the fermion doubling problem and allows to intro- 
duce suitable boundary conditions at the system edges. 
To illustrate its utility we consider boundary conditions 
mimicking the zigzag edges of a graphene sheet .'^ In- 
spired on the procedure proposed by Tworzydlo and 
collaborators,'^^ we develop a method to compute trans- 
port properties of massless Dirac particles in systems 
with different edge boundary conditions. In addition, 
the method we put forward also allows to account for 
the presence of an external magnetic field. 

We apply the method to compute the average conduc- 
tivity at the charge neutrality point of large disordered 
graphene systems. We consider diagonal disorder and 
compute the average conductivity by taking disorder en- 
semble averages. We scale the conductivity with the sys- 
tem size to numerically obtain the corresponding univer- 
sal curve, a matter of an intense controversy in the recent 
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literatureP'i^l'l^ We also study the conductance fluctua- 
tions at the charge neutrality point, providing numerical 
evidence of a non universal behavior. 

The presentation is organized as follows. In Sec. |TT] 
we review the discretization strategies used to solve the 
Dirac Hamiltonian. In Sec. Illll we construct a discretiza- 
tion procedure to calculate the Landauer conductance 
of two-dimensional massless Dirac particles confined to a 
ribbon with boundary conditions that mimic zigzag edges 
in graphene. Ballistic transport is discussed in Sec. |IV[ 
where we apply the method to reproduce the standard 
results found in the literature. In Sec. [V] we apply the 
transfer matrix method to study diffusive transport in 
long-range disordered graphene stripes. Next, we show 
how to include a magnetic field in the method and com- 
pute the conductance at the (anomalous) quantum Hall 
regime in Sec. |VI[ F inally, we present our conclusions and 
outlook in SecJVlIl 

II. MASSLESS DIRAC EQUATION IN A 
LATTICE 

Let us consider the two-dimensional massless Dirac 
Hamiltonian 

H = -ihv{(T^dx + aydy) + U{r), (1) 

where v is the velocity of the massless Dirac fermions, 
ax and ay are Pauli matrices, and U{r) is the potential 
energy landscape. The eigenvalue problem reads 

H-^ = , (2) 

where E is the fermion energy and is a two-component 
(spinor) wave function of the form 

(3) 

For graphene, the two components of the spinor ^' cor- 
respond to the electron wave function amplitudes at each 
of the two intercalated triangular lattices, denoted by A 
and B, that form the honeycomb lattice. Equation ([ij is 
the effective low energy graphene Hamiltonian for elec- 
trons in the vicinity of the Brillouin zone -fiT-point, called 
if- valley. By replacing ay — > —ay one obtains the corre- 
sponding equation for the K'-val\ey^ 

In the case of topological insulators, the spinor com- 
ponents of vjf are interpreted as the electron z-axis spin 
projection. 

Our aim is to find a suitable lattice representation for 
the H operator in Eq. ([T]). The square lattice, that is 
very efhcient to discretize a two-dimensional Schrodinger 
Hamiltonian, would be the most natural candidate. How- 
ever, it is well established in the lattice gauge theory 
literatur^i^Ilillll ^jj^t this choice leads to wrong results. 

Let us understand this statement and motivate the 
discretization scheme we propose by discussing the typi- 
cal problems that arise in the discretization of the Dirac 



equation. For the sake of simplicity, let us consider the 
simplest possible setting, namely, a one-dimensional sta- 
tionary massless fermion field. The Dirac equation reads 

- ihva^dx'i'ix) = E'${x) . (4) 

Despite the simplicity of Eq. Q, its solution by dis- 
cretization requires some care. The naive evaluation of 
the derivative as = -I- A) - *(a::)]/A is un- 

suited, since it produces a non-Hermitian Hamiltonian. 
The reason is that the left and the right hand sides of 
Eq. Q are evaluated at different points, namely x and 
x + A/2. The evaluation of d^^ix) in Eq. Q by tak- 
ing dx'^ix) = [^{x + A) - *(x - A)]/2A also fails. In 
this case, one finds an Hermitian Hamiltonian, but with 
a wrong spectrum: The number of fermions in the first 
Brillouin zone is doubled. This is the so-called fermion 
doubling problem. It arises because the lattice size A is 
half of the spacing between the points used to estimate 
dx^ix)^ 

To correctly solve Eq. Q by discretization, it is nec- 
essary to evaluate dx'^{x) and '^{x) at the same mesh 
position and to use a single effective lattice spacing. Be- 
low, we review two different strategies that satisfy these 
requirements. 

Stacey discretization^^ The fermionic fields are deter- 
mined at the lattice points x^ — m/S. and the Dirac 
equation is evaluated at the lattice midpoints, namely, 
mA-f A/2. 

In other words, since the derivative is dx'^{x + /S./2) = 
-I- A) - vI'(a;)]/A and both sides of Eq. Q must 
be evaluated at the same position, the fermionic fields 
are approximated by the average between adjacent sites, 
namely, -I- A/2) = [^{x + A) + ^{x)]/2. 

The resulting system of coupled finite difference equa- 
tions reads 

_ / V^m+l - V^m ^ +t"'\ (5) 

A V Vm+l - V'm / 2 V Ipm+l + V'm / ' 

where ipm = V'('7^A) and VVt = ipimA). 

Susskind discretization!^ This scheme is based on a 
slightly more involved construction than the one just pre- 
sented. The fermionic fields, ip and ip, are defined in two 
intercalated lattices. The lattices have a lattice parame- 
ter A and are offset by a distance A/2. 

Let us consider the spinor component ip to be de- 
fined at Xm = rnA and the component ip defined at 
Xfn — rnA — A/2. Given that a^ is non diagonal, one 
can evaluate the upper component of Eq. Q at x and 
the lower component at a; — A/2. 

The Dirac equation, Eq. Q, is approximated by the 
coupled finite difference equations 

A V - V'm-^l / V ^rn J ' 

where V'm = ip{mA) and ipm = tp{mA — A/2). No- 
tice that the definitions of ipm and ipm for the Susskind 
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discretization differ from tlie ones used in the Stacey dis- 
cretization scheme. 

In summary, both reviewed discretization procedures 
fulfill the desired goals: They render Hcrmitian Hamilto- 
nians with the correct spectrum at the continuum limit, 
eliminating the fermion doubling problem. 



III. TRANSFER MATRIX DISCRETIZATION 

In this section we describe how to combine the trans- 
fer matrix method and the discretization schemes pre- 
sented above to compute the Landauer conductance for 
graphene ribbons with zigzag edges. We use the Susskind 
discretization along the ribbon transversal direction, be- 
cause it allows for the implementation of zigzag-like 
boundary conditions at the ribbon edges. Along the lon- 
gitudinal direction, we use the Stacey scheme, since it can 
be nicely combined with the transfer matrix method. 

We consider a geometry where source and drain are 
placed at two opposite sides of a rectangular sheet. The 
graphene sheet has a length L in the direction connecting 
source and drain, that is longitudinal to the electronic 
propagation, and a transversal width W . 

We take the x-axis along the longitudinal direction, 
and write Eq. ([I]) as 



(7) 



where V ~ {U — E)A/hv. The Landauer conductance is 
obtained by solving Eq. ([t]) in a lattice using the transfer 
matrix method. For notational convenience, the lattice 
size A is accounted for in the definition of V, a dimen- 
sionless quantity. 

The graphene zigzag-like boundary conditions are im- 
plemented by making ip = ai y = and ■0 = at 
y — W, corresponding to opposite edges of the ribbon.'^i' 
In Appendix]^ we show how to implement the finite dif- 
ference method to address armchair-like boundary con- 
ditions. 

The discretization of Eq. ([t]) requires the use of, at 
least, 3 intercalated lattices. The discretization we put 
forward in this study is shown in Fig. [T] We recall that, 
for graphene, the spinor components and ^ stand for 
(continuum) wave function amplitudes at the sublattices 
A and B, respectively. Accordingly, the discrete repre- 
sentation of the spinor component tp is given by the value 
of ■0 at sites of the lattice A (indicated by filled circles) . 
Likewise, tp is discretized by taking its value at the sites of 
the lattice B (open circles) . The Dirac equation is evalu- 
ated at the auxiliary intercalated lattice, represented by 
crosses in Fig. [ip^l 

The finite difference representation of the partial 
derivatives appearing in Eq. ^ is given by 
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FIG. 1: Representation of the discretization lattice. The tp 
spinor is evaluated at the positions corresponding to the filled 
circles, while tp is evaluated at the open ones. The crosses 
represent the intercalated lattice where the Dirac equation is 
evaluated. 



Here, d^ip and dytp are evaluated at the site s of the auxil- 
iary intercalated lattice, located at the center of a square 
cell whose corners correspond to the lattice A sites where 
i/j is evaluated. The diagonal arrows in Fig. [l] illustrate 
the input required to compute the partial derivatives of 
-0 at a given site s (indicated by a cross) of the auxiliary 
lattice. 

Analogously, the expressions for the lattice B read 

dxi^ =^(V'm-l-l,n - ■0m,«) + (V'm+l.n-l " V'm.n-l) (9) 



■0m+l,n-l) + ilpm,n ~ '0m,n-l) ■ 

Finally, the lattice approximation for V"^ in Eq. ^ is 



VtP^V^, 



(10) 



which is evaluated, as above, at the site s of the auxiliary 
lattice. More precisely, Vm.n gives the value of V at the 
midpoint between the sites (m, n) and (m + l,n) of the 
lattice A. The evaluation of Vtp is represented by the 
horizontal arrows in Fig. [T] 
For the spinor component tp 



0, 



m.n 



(11) 



gives the value of Vtp at the auxiliary lattice site s at the 
midpoint between the sites {m,n) and (m -I- l,n) of the 
lattice B. 

Throughout this paper, we adopt the convention 
that the m index is associated with the a:-axis, while 
n labels the mesh points y-coordinate. For nota- 
tional convenience, we define the column vectors xp,n = 



("^m,!,--- ,i>m.NY' and Xpra = (V"™,!, 

also introduce the column vector 



^ r 



tpn 



,ipm,N) ■ We 



(12) 
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that represents the spinor of Eq. ([s]) evaluated at both 
A and B sublattices for a given m, corresponding to a 



X = mA. In Fig. [T] single columns represent spinors 
of different m's. 



Plugging ([8| to ([TT]) in Eq. ([?]) we obtain 

(I + N+){^P„^+l - i>rn) =i{i - N+){^^>„^+l + ^m) " + (13) 

(I + N")('0„i+i - ■0„j) =i(I - N")(i/J„j+i + ■0„i) - i¥,„('0,„+i + ■0„j) 

where I is the N x N identity matrix, and are super- and subdiagonal matrices, namely, (N"'")„„' = Sn,n'~-i, 
(N~)„„' = (5„,„'+i, while V™ and V™ have matrix elements given by (¥„)„„' = VmnSnn' and (¥„)„„' = VmnSnn' ■ 



By writing Eq. 13 in terms of and '4'm+i, 

N+) iV" 



(I + N+) 



(I + N-)-i(I-N-) 



(I + N-) + i(I 



r 



(14) 



allows us to identify the transfer matrix M defined by 

= M„*„, (15) 

namely, 

I - iX„,, 



where 



and 



-(I-N+) V" 



-(I -N- 



: + N+ 

1 + N- 



(16) 



(17) 



(18) 



We define the current operator in the a;-direction as 



1 

= r 



J_ 
J+ 



(19) 



with J+ = I + N+ and J_ = I + In Appendix |b] 
we show that this choice of Jj; guarantees charge conser- 
vation. As a consequence of the zigzag boundary condi- 
tions, our discretization method does not preserve sym- 
plectic symmetry. This is in contrast with the method 
put forward in Ref. |151 that deals with periodic bound- 
ary conditions along the y-direction and treats infinitely 
wide systems. 



A. Landauer conductance 

The transfer matrix of the whole graphene sheet can, in 
principle, be computed using M = nm=i ^m- However, 
it is well known that this scheme is numerically unstable: 
It produces both exponentially growing and exponen- 
tially decaying eigenvalues. Limited numerical accuracy 



prevents one from retaining both sets of eigenvalues .^^^ 
The numerical instability problem is solved by working 
with compositions of unitary scattering matrices ,'^^'^11 at 
the expense of introducing additional diagonalizations 
that slow down the computation. We adopt the scheme 
proposed in Ref. [15] with a small modification, which we 
briefly describe in what follows. 

In order to translate the transfer matrix M into a uni- 
tary scattering matrix 



r t' 
t r' 



one needs to solve the equation 

flight ^ iv[[$'°^* 



with 



right 



tah<l>, 



+ 
b ' 



left 



(20) 



(21) 



(22) 



for each asymptotic propagating channel a at the source 
(left) or drain (right) contacts. We consider the situation 
where all channels are open. This is similar to the set- 
ting where the contacts have a high doping and a large 
density of states to realistically model standard metallic 
contacts,!^ minimizing the effect of a contact resistance. 

In the large wave vector limit, the states (jy^ (a = 
1,2,..., N) are given by the 2A^ eigenstates of the current 
operator Jj. , normalized such that each mode carries the 
same current. That is 



6± = 7-1/2 

a J a 



where and ja are given by 



(23) 



(24) 



with (V'a)n = (i^a)N+l-n- 
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To obtain a closed-form expression for S in terms of 
M, we write M in the basis of channels 0j through the 
similarity transformation 

M„=MM„jM"i (25) 

where 

{R~%Nx2N = {(l^tr- - ,<P%,(l)i,--- .(Pn)- (26) 

This M is not the same as the one used in Ref . [TSl We use 
a different transversal discretization scheme and K is de- 
fined to facilitate the projection of single-contributions to 
the transmission from different channels. The later char- 
acteristic is very useful for the method implementation, 
as justified in the next section. 

The spinor degrees of freedom of M arc separated into 
four N X N blocks, 
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FIG. 2: Band structure of Dirac fermions in a infinite ribbon 
of width W — lOOA with zigzag-like boundary conditions 
obtained using the method described in the text. The energy 
scale is set by fixing A = lOao, where ao = 1.42 A is the 
carbon-carbon bond length in graphene. 



and follow the same strategy employed in Sec. Ill.b of 
Ref. [TSj to write the S matrix in terms of the M sub- 
blocks. Finally, the conductance is evaluated as usual by 
the Landauer formula 

G = ^£^T (28) 

where gs and stand for the spin and valley degeneracies 
and T is the transmission, namely, T — ti[tH]. 

B. Spurious mode and its elimination 

Let us discuss the properties of the transfer matrix for 
the simple case of a pristine graphene ribbon with zigzag 
edges. Due to translation invariance, all transversal slices 
m are equivalent, and the transfer matrix Mpuc = M^, 
corresponding to the primitive unit cell, is sufficient to 
describe the transport. 

For a given energy E, the transfer matrix Mpuc has 
A'ch eigenvalues of the form e**^""^ with ka real. Hence, 
A'ch gives the number of propagating channels and ka is 
associated with the wave number along the propagation 
direction of the corresponding channel. This allows us to 
obtain the dispersion relation for a graphene ribbon. For 
a given ribbon width VF, we vary the energy E in small 
steps and calculate the eigenvalues and eigenchannels of 
Mpuc- The result for the discretization of the X- valley 
is shown in Fig. [2] 

The band structure shown in Fig. [2] coincides with 
the results obtained by the continuum limit treatment 
of the zigzag boundary conditions of Ref. 11 around the 
K valley. As expected, the low energy states are also in 
good agreement with those of the nearest-neighbor tight- 
binding model. The only discrepancy is the appearance 
of an "extra" mode around kA/n sa 1. To understand 



its nature, it is convenient to analyze the corresponding 
eigenstate. This is done next. 

The transfer matrix Mpuc has two zero-energy eigen- 
modes, whose typical transversal wave function ampli- 
tudes are shown in Fig. [3j The "smooth" mode corre- 
sponding to Fig. is similar to the analytical solution 
presented in Ref. JT, In addition, we find a spurious 
mode that oscillates on the scale of the lattice size dis- 
cretization, shown in Fig. [3]d. 

In finite difference calculations, states showing oscil- 
lations at the scale of lattice size are expected to be 
badly described by the discretized equation. Usually, 
they correspond to the states of higher energy in the 
spectrum. The general rule is that the first half or so of 
the computed states are well described by the discretiza- 
tion. Hence, provided one is interested in the low en- 
ergy part of the spectrum, the inaccurate high frequency 
states do not cause any problem. These observations 
do not apply here. For example, in the Tworzydlo and 
collaborators^^, that employ a Stacey-like discretization 
to handle the fermion doubling problem, the effective lat- 
tice size is doubled to evaluate the derivatives of As 
a consequence, the high energy states at the top of the 
spectrum are poorly described and spurious states arise. 
To eliminate these spurious states, very long filters at 
each side of his graphene stripe are introduced. 

In our case, the Susskind discretization preserves the 
lattice size, however the peculiar nature of the edge state, 
characteristic of ribbons with zigzag edges, pulls down 
the energy of the maximally transversal oscillating state. 
The combined oscillations in both sublattices render a 
dispersionless mode in the finite difference representation 
of the momentum operator. 

To eliminate (up to a very good accuracy) the spuri- 
ous state illustrated in Fig. [3]3, we simply eliminate the 
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FIG. 3: (Color online) Spinor ^ projection amplitudes, ?/) and 
V^, corresponding to the zero-energy dispersionless modes of 
Mm: (a) Physical mode and (b) spurious solution. 

corresponding channel in definition of M: 

(M-')2A^x2Ar-2 = (0^, • • • , ■ ■ • , 0^V-l)- (29) 

The inversion of the is taken by transposing the 
matrix and by inverting the normalization of the eigen- 
vectors . 

IV. BALLISTIC TRANSPORT 

In this Section we show that our method reproduces 
known results for ballistic transport of Dirac fermions 
in graphene ribbons with zigzag edges. We consider 
graphene stripes of length L and width W . The potential 
energy is C/(r) = C/q for < x < L. Throughout this sec- 
tion A = IOoq, where — 1.42 A is the carbon-carbon 
bond length in graphene. 

The conductance is computed using the procedure de- 
scribed in Sec. |III A} Figure |4] shows the transmission co- 
efficient T = Tx+Tj^i as a function of the energy E prior 
to the spurious mode elimination. We obtain Tk = Tk' ■ 



The conductance around the charge neutrality point 
E — Uq = is twice the expected value. Instead of a 
single open mode, corresponding to edge state in the 
zigzag boundary,^ we observe two open modes at low 
energies. The "extra" mode is just the spurious one dis- 
cussed above. 

The spurious mode is removed following the recipe in- 
dicated in Sec. [Ill B1 Figure|5]shows that the transmission 
coefficients for the K and K'-valleys are no longer degen- 
erate. At low energies the transfer matrix method gives a 
perfect unit conductance for E > around K (the same 
for E < around K' .) On the other hand, subtract- 
ing the corresponding unit step from Tk and Tk' both 
transmission coefficients coincide and show particle-hole 
symmetry. The total conductance is given by the sum of 
the contributions from both valleys. 

The low energy features of the conductance can be un- 
derstood by inspecting the band structure of a zigzag 
nanoribbon, as shown in Fig. [2] For energies E close to 
the charge neutrality point, after eliminating the spuri- 
ous mode, there is a single allowed energy band. For 
-E > 0, right moving states belong to the if- valley. Con- 
versely, left moving states belong to the isT'-valley. The 
-fC'-valley dispersion relation can be obtained by making 
fc — > — fc in Fig. [2] In our geometry, the drain is at the 
left, while the source reservoir is connected to the right 
contact. Hence, ioi E > in the vicinity of ii^ = there 
is a single open channel, which belong to the if -valley. 
The picture is reversed for E < 0. 

As \E — Uo\ is increased. Fig. |5] shows a sequence of 
steps decorated by Fabry-Perot-like oscillations. The 
dashed lines indicate the energy threshold for successive 
subband openings inferred from Fig. [2j We observe a 
good agreement between the dashed line energies and the 
steps in Fig. [5] The mismatch between the propagation 
energy in the graphene sheet and the contacts gives raise 
to the Fabry-Perot-like interferences. This assertion can 
be quantitatively tested by calculating T for ribbons with 
a fixed width W for different lengths L. We observe that, 
as L is increased, the Fabry-Perot interference oscillate 
more rapidly, scaling with (figure not shown). 
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FIG. 5: Transmission T for graphene ribbons of width W = 
100 A and length L = 300 A as a function oi E - Uo in eV 
using a massless Dirac equation centered at the (a) A'-point 
and (b) /("'-point. The vertical dashed lines are placed at the 
threshold energies of the propagating modes in the depicted 
energy interval. 



Due to their rough edges, such Fabry-Perot interfer- 
ence patterns are unlikely to be observed in lithographic 
graphene nanoribbons. The situation is more promising 
in smooth edge graphene ribbons, made by unzipping 
carbon nanotubes.'^S] 

We note that alternative discretization schemes can 
be used to construct a transfer matrix to compute the 
transport of Dirac fermions in ribbons with armchair- 
like boundary conditions,!^ and possibly also for chiral 
boundary conditions.'^ Our motivation to focus only at 
the zigzag case is twofold: (i) In distinction to the others, 
this boundary condition does not mix the valleys isospins. 
This allows to solve the problem in the K and /^''-valleys 
separately, which speeds up the computation by a factor 
about 2^. {ii) The method was put forward to compute 
transport properties of realistically large samples in the 
presence of long range disorder, for which the sample 
boundaries should not play a key role. 



V. MESOSCOPIC EFFECTS 

The electronic transport in graphene at the charge 
neutrality point has been the subject of intense exper- 
imental and theoretical work, as reviewed, for instance 
in Refs. I4|5I While it is by now established that in good 
mobility graphene samples the observed conductivity a 
is of the order of e^//i, its value still depends on sample 



disorder. In the diffusive regime, where the system size 
L is much larger than the electron mean free path £, a 
is expected to show a weak system size dependence. In 
the ballistic case, where L < £, the conductivity depends 
strongly on the system geometry and it is more appropri- 
ate to describe transport properties by the conductivity. 
It is likely that the graphene samples used in many of the 
reported experiments correspond to a ballistic-diffusive 
crossover .SI 

In this Section we apply the finite difference method 
to compute the graphene conductivity a at the charge 
neutrality point in the presence of long range disor- 
der. We obtain the expected anti-weak localization re- 
lation between a and the system size,'^^ namely, (a) = 
{gsgvC^ /TTh)\n{L/£). Next, we study the conductance 
fluctuations at the charge neutrality point and show that 
they are not universal. 

The transfer matrix approach gives the two-terminal 
conductance G. The conductivity is obtained from the 
relation a — {L/W){G), where (• • • ) stands for a disorder 
ensemble average. Since we assume that there is no inter- 
valley scattering, there is no need to choose a finite corre- 
lation length for the potential fluctuations,^''' in distinc- 
tion to other previous theoretical analysis.* 6 7 17 is 26 28] 
As in Ref. [15] we let the potential at each lattice point 
fluctuate independently. More specifically, we draw 
Um,n from an uniform distribution within the interval 
{-5U,SU]. 

We compute the conductivity for 10'^ realizations for 
different values of 6U^ corresponding to different mean 
free paths. We assume that the conductivity follows the 
relatiorfl^ 



— = - In 

Go TT 



L 



C{5U) 



(30) 



where Go — gsgv^^/h. Here the first term at the right- 
hand side is just the weak anti-localization relation with 
1*{SU) associated with £, while the second one accounts 
for finite size corrections and/or a contact resistance. For 
every set of disorder realizations with a given value of 6U 
we vary the system size {L = W — 30, 60, 90, and 120A) 
and compute {G[L)). With the help of (30), these data 
are used to find t{5U) and G{5U). 
Figure [6] shows 



In 



L 



(31) 



obtained from the simulations described above. We re- 
mark that our method allows to explore a much larger 
range of L/f than any previous study. 

The current understanding of the conductance fluc- 
tuations is less clear. Analytical works make predic- 
tions of universal conductance fluctuations (UCF) in 
graphene.'22l3SI However, since these results rely on di- 
agrammatic expansions in powers of {kFt)~^, they are 
applicable for large doping but are hardly relevant to 
explain the fluctuations at the vicinity of the charge neu- 
trality point. The available experimental datsPS do not 
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FIG. 6: Conductivity a at the charge neutrality point as a 
function of system size L/C. Points correspond to differ- 
ent system sizes L = W = 30, 60, 90, and 120A and disorder 
strengths SU = 2, 4, 6, and 8. The dashed line is given by 
Eq. (3i|. The lattice spacing is A = 5ao. 



provide a very consistent picture, possibly indicating an 
absence of UCF in graphene at zero doping. 

An early numerically worlP studied conductance fluc- 
tuation in graphene in the presence of short range dis- 
order. More recently, a systematic analysis of the 
conductance fluctuations in graphene sheets has been 
performed.'^ The study reports conductance fluctuations 
as a function of doping, system size L/^, and disor- 
der strength. The authors consider graphene sheets 
with transversal periodic boundary conditions (no edges) 
and use the finite difference method put forward by 
Tworzydlo and collaborators.^^ 

We analyze the variance of the conductance fluctua- 
tions, var(G) — (G^) — (G)^, at the charge neutrality 
point as a function of L£* . In distinction to Ref. UHl we 
find no evidence of universal conductance fluctuations, 
even going deep into the diffusive regime L/i* ^ 1. As 
stated above, we believe that the strong discrepancy be- 
tween our results a nd the prediction from diagrammatic 
perturbation theorypSEHl jg (j^e to the breakdown of the 
perturbation theory when kp£ 1, that is, when one ap- 
proaches the charge neutrality point. In such situation, 
numerical simulation such as the one presented here give 
the correct way to approach the problem. 



VI. INCLUSION OF A MAGNETIC FIELD 

In this section we show how to modify the transfer ma- 
trix discretization scheme presented in Sec. |III| to account 
for the presence of an external magnetic field. This allows 
one to calculate (in the single-particle picture) the con- 
ductance of ballistic and long range disordered graphene 
sheets for arbitrary values of a perpendicularly applied 



FIG. 7: Conductance fluctuations variance var(G/Go) as a 
function of system size L/£*. 



magnetic field. 

As standard, an external magnetic field is introduced 
in the Dirac Hamiltonian by minimal coupling.I^ In what 
follows, we consider the vector potential A ~ {—By, 0, 0) 
that corresponds to the constant magnetic field B — Bz. 
It is important to mention that, since the discretiza- 
tion schemes along the x and y directions are different, 
gauge invariance is broken and the choice of A has to 
be cautiously made. The minimal coupling substitution 
p ^ p {e/c)A modifies Eq. Q to 



{d^ - ieBy/hc)-^ = -~ii<Jzdy + a^V/A)'^ 



(32) 



Note that the magnetic field does not couple different 
valley pseudospin projections. Equation (32 1 refers to 
the A'-valley. As before, by replacing ctj, — > —ay one 
obtains the corresponding equation for the i^'-valley. 
The transfer matrix discretization for Eq. ( 32 ) is con- 



structed in the same way as the one presented in Sec |III| 
The expression for is given by Eq. ( 16 1, but now X„ 
is given by 



2<po 



200 



(33) 

where (j) = BA^ is the magnetic flux through a square 
unit cell (lattice A or B), 0o = h/2ec is the magnetic flux 
quantum, while D and D are diagonal matrices, whose 



nS„ 



and 



matrix elements read D„„ 
l/2)^„„'. 

To illustrate the results that can be obtained by this 
method, we compute the Landauer conductance of a bal- 
listic graphene sheet in the anomalous quantum Hall 
regime. 

Figure |8] shows the if-valley conductance for a rect- 
angular graphene ribbon geometry, where L — 400 A, 
W — 200A and A = oq, subjected to different perpen- 
dicular B fields. The energy axis is scaled by {Bq/BY^^, 
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where Bq = 20 T to make evident that the Landau 
level energies scale with as expected for Dirac 

fermions.'i' More precisely, Ej^ — hucN^^"^, where the cy- 
clotron frequency ujc = vp(2eB/hY^^ . The transmission 
steps in Fig. [8] nicely coincide with the Landau level en- 
ergies. The different B-field conductance curves strongly 
overlap. We observe that as E is increased, the conduc- 
tance steps become less well-defined. This behavior is 
particularly evident for J3 = 20 T and become less pro- 
nounced as B is increased. The interpretation is that 
with increasing B, the edge states become more local- 
ized and hence tunneling between edges, which destroys 
quantization, is suppressed. 




0.02 0.04 0.06 0.08 0.1 



{E-Uo)/{B/Boy'^ [eV] 

FIG. 8: Transmission T for different magnetic field values as 
a function of E ior L = 400ao and W — 200ao. To make the 
magnetic field dependence of the Landau levels in graphene 
expficit, the energy axis is scaled by (Bq/ B)^^^ , with Bo = 
20T. After scafing all curves largely coincide. For the sake of 
clarity, the transmission for B = 40 T is offset by +2 units, 
while T for B = 80 T is offset by +4 units. 



graphene and surface states in 3D topological insulators. 

The use of a finite difference description of the Dirac 
equation to calculate transport coefficients under given 
boundary conditions demands tailormade discretization 
schemes. We consider a system with a rectangular ge- 
ometry. In the longitudinal direction, parallel to the axis 
connecting source and drain, a Stacey-like discretization 
allows us to implement the transfer matrix formalism. In 
the transversal direction, a Susskind-like discretization 
makes possible the inclusion of zigzag-like edge bound- 
ary conditions. By construction, our method avoids the 
fermion doubling problem. 

The proposed method becomes numerically advanta- 
geous with respect to the standard ones, that employ an 
atomic basis, when treating systems where long range 
disorder is dominant. In this case, the disorder range f 
sets the discretization length scale. That is, A < ^ in- 
stead of A ~ flo as in any atomistic basis discretization. 
For a ribbon of length L and width W, the computation 
time of transport properties scales with {L/A) x {W/ A)^. 
Hence, for a given geometry, our discretization method is 
of the order of (^/ao)^ times computationally faster than 
standard atomistic ones. 

To illustrate how the method works, we analyzed the 
transport properties of ballistic and diffusive graphene 
ribbons in the absence of magnetic field, as well as the 
anomalous quantum Hall plateaus in a graphene sample 
under a strong perpendicular magnetic field. 

We believe that the most promising features of our 
method is that it allows considering different kinds of 
boundary conditions. The present study addresses the 
case of zigzag-like edges, but the transversal discretiza- 
tion can be modified to mimic arm-chair (see Appendix 
|A| ) and chiral edges, as well as non-trivial periodic bound- 
ary conditions suitable to describe surface states of cer- 
tain topological insulators. We believe our method can 
be very helpful for the study of issues related to the in- 
triguing nature of boundaries in Dirac systems. 



For the system we analyze, the K'-valley transmis- 
sion coefficient ioi E - Uq > is Tk' ^ Tk - I- The 
corresponding figure is not shown, but bears a strong 



analogy to the situation we discussed in Sec. IV The 
full conductance is given by the contribution of both 
K and K' valley pseudospin components, namely G — 
{gse^ /h){TK -t-T^i) and consistent with the result found 
in the literature.'^J^ 
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Appendix A: Armchair graphene ribbons 



VII. CONCLUSION AND OUTLOOK 

In this paper we put forward a numerical method to 
calculate transport properties of two-dimensional mass- 
less Dirac particles. The numerical procedure combines 
a finite difference description of the Dirac equation with 
a current conserving transfer matrix formalism. The 
method is useful to study the Landauer conductance in 



In this Appendix we modify the discretization proce- 
dure presented in Sec. |III| to address graphene ribbons 
with armchair-like edges. To satisfy the boundary condi- 
tions of such systems on e ne eds to mix both K and K' 
Dirac- valley components,!^^ namely. 
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The armchair-like boundary conditions at the graphene 
ribbon edges are accounted for by taking^ ^ 



*(r) 



x=0 



*'(r) 



x=0 



and 



*(r) 



iSKW 



x=W 



*'(r) 



x=W 



(A2) 



(A3) 



where 5K — 47r/(3ao) and oq is the graphene lattice con- 
stant. Here we chose the x-axis as the transversal direc- 
tion and the y-axis as the electron propagation direction. 
Equations ( A2 ) and ( A3 ) guarantee that the wave func- 



tion amplitude vanishes on both sublattices at the ribbon 
edges X = and x = 

In distinction to the zigzag case, described in terms 
of a single- valley Dirac equation given by Eq. ([T]), here 
we have to solve the effective 4-component spinor Dirac 
Hamiltonian 











<7xdx 



aydy 



= E 



(A4) 



with boundary conditions given by Eqs. (A2) and (A3 1. 

The problem is cast in the discrete space as follows: 
The fields ^E" and ^P' are evaluated at 



1 





= nA) 


n 


= 0, • 


• .N 




= nA) 


n 


= 1,- 


■ ,N 




= nA) 


n 


= 1,- 


■ ,N 




= nA) 


n 


= 0,- 


■ ,N 



1 



(A5) 



where A is the lattice spacing and m labels the mesh 
position, y = to A. The armchair-like edge boundary 
conditions ( A2 ) and ( A3 1 read 



i^mix = iVA) = e^^^'^CGx = iVA), 
ip^ix = 0) = ^^„(a; = 0), 

^Ux = NA) = e-™Vi„,(a; = iVA), 



(A6) 



where W = (iV+ 1/2)A, as in Ref. El We set A = oq to 
make easy the comparison with the tight-binding model 
results.'^ 

Following the steps of the discretization procedure de- 
scribed in Sec. Ill Eqs. ( A5 ) and ( A6 ) lead to a modified 



-matrix given by 




-{I-'N- 






B 







(A7) 

where A and B are matrices whose single non-zero ele- 
ments are Ajvjv = —e~'^^^^ andBi^i = —1, respectively. 



Similarly, the J-matrix reads 



/I + N+ 







+ N- 
At 




A 





(A8) 



The transfer matrix M is defined, as standard, by 



(A9) 



As in Sec. |III) is related to X„i by Eq. (|16|. 

Due to translational invariance, the transfer matrix of 
pristine armchair graphene ribbons does not depend on 
the index to, that is Mpuc = M^. For a given energy 
E, Mpuc has a subset of eigenvalues of the form e''"'"^, 
with ka real. Those correspond to propagating modes 
and allow one to infer the system band structure. For a 
detailed discussion see Sec. IIIIBI 
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FIG. 9: Dispersion relation E/t as a function of kaof-K for 
graphene ribbons with armchair- like edges, where t is the 
graphene tight-binding nearest-neighbor hopping integral, (a) 
Insulator ribbon, TV = 24. (b) Metallic ribbon, = 25. 

The described discretization scheme reproduces with 
good acuracy the low-energy nearest neighbor tight- 
binding dispersion relation for graphene ribbons with 
armchair edges. As N is increased, we obtain the cor- 
rect sequence of two insulator band structures followed 
by a single metallic one. Figures and [9]d show the 
band structure for typical metallic and insulator ribbons. 
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respectively. As in the zigzag case, there are spurious 
modes in the metaUic ribbon, whose wave functions os- 
cihate on the scale of the lattice spacing. Those must be 
eliminated to understand the system transport proper- 
ties. This can be done in a similar way as that discussed 
in Sec. HITbI 

Figure [TO] shows the dependence of the lowest confined 
state energy versus ribbon width given in terms of N. For 
comparison, the corresponding energy obtained from the 
nearest-neighbor tight-binding model as also shown. The 
agreement is very good, comparable to the one obtained 
solving the Dirac equation without discretization.!^ 



graphene effective Dirac Hamiltonian, the transfer ma- 
trix method can accurately describe the band structure 
of graphene ribbons with armchair edges. After the elim- 
ination of the spurious modes that appear in the metallic 
case (not done here), the computation of the conductance 
follows the procedure described in Sec. 

Further extensions of the method to deal with other 
kinds of edge boundary conditions, such as the chiral 
ones, is also possible. The later requires a modification of 
both the primitive unit cell and the boundary equations, 
namely, (A5) and (A6). 
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Appendix B: Charge conservation. 
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FIG. 10: Lowest energy E* /t of the graphene ribbons with 
armchair-like edges as a function of the ribbon width W = 
Nao- Crosses correspond to the energies values calculated 
using the finite difference method, while the dots stand for 
nearest-neighbor tight-binding model results. 



In order to conserve the total current at any arbi- 
trary ribbon cross section, the discretized current opera- 
tor must fulfill the following condition 



(*„+l|Jx|*m+l) = (*m|Jx|*.r 

As a consequence, one writes 

Mt„JxM„ = J„ 
which is equivalent to 



(Bl) 
(B2) 

(B3) 



Using Mm as given by Eq. (16), we show that in order 



to preserve current, Ja, must satisfy 

J-ixJ„J, = X,„. (B4) 
Indeed, this can be readily shown with the help of Eqs. 



These results show that, by using the full 4-component (171, (181, and ( 19 1 
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